Cite this Article
Ouyang Zhenyu, Lin Jianzhong, Ku Xiaoke. Dynamics of a self-propelled particle under different driving modes in a channel flow. Chinese Physics B, 2017, 26(1): 014701  
Permissions
Dynamics of a self-propelled particle under different driving modes in a channel flow
Ouyang Zhenyu, Lin Jianzhong †, Ku Xiaoke
State Key Laboratory of Fluid Power and Mechatronic Systems, Zhejiang University, Hangzhou 310027, China

 

† Corresponding author. E-mail: mecjzlin@public.zju.edu.cn

Abstract

In this paper, a model that combines the lattice Boltzmann method with the singularity distribution method is proposed to simulate a self-propelled particle swimming (exhibiting translation and rotation) in a channel flow. The results show that the velocity distribution for a self-propelled particle swimming deviates from a Maxwellian distribution and exhibits high-velocity tails. The influence of an eccentric potential doublet on the translation velocity of the particle is significant. The velocity decay process can be described using a double exponential model form. No large differences in the velocity distribution were observed for different translation Reynolds numbers, rotation Reynolds numbers, or regular intervals.

PACS: 05.20.Dd;47.11.Qr;47.57.E-
Keyword:lattice Boltzmann method;singularity distribution method;self-propelled particle;driving modes
1. Introduction

Self-propelled particle swimming has gained wide attention due to its relevance to biological and technological applications. Micro-swimmers may use various driving mechanisms, such as flagellar propulsion, [1–3] beating cilia, [3] surface distortions, [4] chemical reactions, [5, 6] or actin-tail polymerisation. [7, 8] Although notable differences between these mechanisms exist, the universality in hydrodynamic behaviour is obvious. For many biologically relevant systems, it is typical that self-propulsive swimming contributes to non-equilibrium dynamics systems. [9, 10] The experimental observation that swimming particles enhance hydrodynamic diffusion has been confirmed by numerical simulations. [11, 12] The presence of swimming particles always affects the stability of suspensions. [13, 14] On the other hand, the rheological characteristics of suspensions that result from swimming particles may play a major role in channel flow. Bioconvection can occur when swimmers share a mean upward movement because of the effects of gravitational torque, chemotaxis, or phototaxis. [15] Therefore, investigation of the motion characteristics and effects of self-propelled swimmers on the properties of fluid media is important to understanding the transformation mechanism of biological system energy and improving the design of micro bionic robots.

Most of the research on the motion of self-propelled particles suspended in Newtonian fluids have employed analytical solutions or numerical methods to examine motion under the condition of low-Reynolds-number flow. From a microscopic perspective, Lighthill [16] first proposed the model of a squirmer, which has a spherical shape with surface tangential velocities, based on a solution of the Stokes equation in a spherical coordinate system. Blake [17] later extended this model to the dynamics of ciliary propulsion of a microscopic organism. Magar et al. [18] developed a model for nutrition uptake by a self-propelled steady squirmer. Recently, Ishikawa et al. [19] developed a database for an interacting pair of squirmers from which one can easily predict the motion of self-propulsion suspensions using a boundary element method. They used the interaction model to study the rheology of a semi-dilute suspension of bottom-heavy and non-bottom-heavy squirmers [20] and the diffusion of swimming micro-organisms in a semi-dilute suspension. [21] They also studied the orientation order in a concentrated suspension of spherical microswimmers. [22] This squirmer model is an accurate combination of an analytical solution and a boundary element method which is precise in dealing with a complex boundary and near-wall flow. However, this method requires a number of integrals and even a singular integral, which requires huge computing resources. Taking into consideration the interaction between particles, Hernández-Ortiz et al. proposed a model in which a force dipole is drawn into a fluid to describe the effects of self-propelled particles. [13, 23] Saintillan and Shelley [24] proposed a model based on the slender-body theory in which swimming particles are regarded as self-propelled rigid slender rods that impose a given shear stress on the surrounding fluid. These models seek to describe the movement mechanism of self-propelled particles and how they affect the characteristics of the surrounding fluid.

The lattice Boltzmann method (LBM), which is widely used in the simulation of various types of flow, [25–29] relies on the fact that the macroscopic dynamics of fluid flow does not depend on the microscopic details of the fluid field. The LBM was first proposed by Ladd et al. to simulate passive suspensions at finite Reynolds numbers. [30–32] Aidun et al. [33–35] developed this method and added ‘virtual nodes’ to the solid boundaries to simulate passive particles near contact. Recently, Ramachandran et al. [36] proposed an LBM model for active suspensions, in which two types of particles (movers and shakers) are defined by different distributions for point forces as a force dipole and a D2Q7 LBM model is used to simulate individual and multiple self-propelled colloidal particles. The researchers found that the velocity of the mover was directly proportional to the asymmetry of the position of the force distribution and that the statistical properties of the velocity field generated by the activity deviated from those of equilibrium fluids. Alarcón and Pagonabarraga [37] proposed a three-dimensional LBM for suspensions of self-propelled particles using a squirmer model in which the appropriate boundary conditions for the Stokes equation at the surfaces of spherical particles were imposed to induce a slip velocity between the fluid and the particles. They employed this method to characterise the hydrodynamic stresses induced by active swimmers, and they found that active stresses played a dominant role in decorrelating the collective motion of squirmers and that contractile squirmers formed significant aggregates.

Several numerical studies have been performed in the past few years to explore the mechanism of self-propulsive particle movement using the lattice Boltzmann method, in which a force dipole or a slip velocity is exerted on the boundaries of the particles as the driving mode. However, it is difficult to apply this approach to particles with irregular geometric shapes. Chwang and Wu [38–40] proposed a singularity distribution method, based on the theory for low-Reynolds-number flows, with which the real movement characteristics of self-propelled particles with complicated boundaries can be simulated. The Stokes flow around an object with arbitrary shapes can be simulated by properly adjusting the type, position, and intensity of fundamental Stokes flow singularities in the interior of an object using the singularity distribution method. In this paper, we present a model that combines the lattice Boltzmann method with the singularity distribution method to simulate a self-propulsive particle swimming in a channel flow. The behaviours of self-propelled particles’ translation, rotation, period of transition between translation and rotation, and decay of the particle velocity were analysed.

2. Numerical methods
2.1. Lattice Boltzmann scheme

The LBM is different from schemes based on the continuum assumption of fluids. In the LBM, the discrete microscopic velocity space is given as follows: [41]

(1)
where f i is the density distribution function; is the streaming velocity in the ith direction in the phase space; i is the collision operator; and the Bhatnagar–Gross–Krook (BGK) collision term is given as follows: [42]
(2)
where is the local equilibrium distribution, τ is the single relaxation time. For a square or cubic lattice, the equilibrium function is as follows: [43]
(3)
where ρ and are the macroscopic fluid density and velocity, respectively:
(4)

In the nine-bit LBGK model, the two-dimensional velocity in the phase space is discretised in nine directions:

The kinematic viscosity for the nine-speed model is and is the lattice speed. In Eq. (3), is equal to 4/9 for i = 1–4 and 1/36 for i = 5–8.

Within the limit of long wavelengths, the lattice Boltzmann equation recovers the following quasi-incompressible Navier–Stokes equation by the Chapman–Enskog multi-scaling expansion: [44]

(5)
(6)
where is the fluid pressure, is the sound speed with and is the dynamic viscosity.

Ladd [31, 32] and Aidun et al. [34] used the momentum exchange method to propose a modified bounce-back rule for a moving wall. If we place the boundary nodes on the links connecting the interior and exterior nodes,

(7)
where denotes the post-collision time; i is the incident direction; is the reflected direction; is the velocity at the particle surface; , where is the translational velocity of the mass centre of the particle; is the angular velocity of the particle; and where is the position of the mass centre. The force and torque exerted by the fluid at are given by the following equations:
(8)
(9)

2.2. Self-propulsion model

For the case in which the inertial effects are negligible in comparison to the viscous forces, under the condition of low-Reynolds-number flow, the Navier–Stokes equations are usually simplified to the Stokes equations

(10a)
(10b)
where is the velocity vector, p is the pressure, μ is the constant viscosity coefficient, is the external force per unit volume, and is the position vector in a three-dimensional Euclidean space. If the value of is zero and the pressure p is constant, the solutions of Eq. (10) are the same as the solutions of classical potential flow in which the velocity vector and the velocity potential satisfy the following harmonic functions:
(11a)
(11b)
We can draw the conclusion that any velocity field of incompressible potential flow with a uniform pressure field ( constant) is the solution of Stokes flow, which means that the two types of flows are similar in kinematics although different in dynamics. Thus, we can add constant to the basic solution of unbounded potential flow as the special solution of Stokes flow that we need. In the spherical coordinate system, the velocity potential that satisfies the Laplace equation can be solved as follows: [45]
(12)
where is the first kind of associated Legendre polynomials and are two arbitrary constants. In the formula, the n = 0 item represents the velocity potential of a point source, and the n = 1 items,
(13a)
(13b)
(13c)
represent the velocity potentials of the doublet for the displacement axis in the , and y directions, respectively. We can obtain the velocities of the point source as follows:
(14a)
(14b)
where i is the direction of the radius vector (the point source is located at the origin of the coordinate system).
(15)
represents the volume flux that flows out from any surface enclosing the origin, which is also called the intensity of the point source ( denotes a point sink). A point source and a point sink can form a potential doublet :
(16a)
(16b)
where D j represents the intensity of the potential doublet for a displacement axis in the j direction. Note that we have not mentioned the other fundamental singular solutions of Stokes flow, including a Stokeslet and its derivatives, known as rotlets, stresslets, and higher-order poles derived from them. The reason is that research has shown that the fluid field driven by self-propelled particles is not like the flow of a Stokes doublet but rather more like that of a potential doublet. [46]

3. Results and discussion
3.1. Computation parameters

In the simulations, we modelled an active self-propelled particle in a typical two-dimensional static flow field, as shown in Fig. 1. The x–y coordinate system was fixed at the centre of the cylinder, with x in the channel wall direction and y in the transverse direction. No-slip and no-penetration boundary conditions were imposed on the channel wall. The particle was 10 in diameter ( . The deformation and bending of the particle caused by the flow were ignored, and the ratio of the particle density to that of the flow was taken to be equal to 1. The lattice spacing was 1 m, and the time step was 1 s. Before we began a simulation, we arranged the periodic boundary length in the flow direction from 10D to 40D to ensure that the effects of the periodic boundary conditions on the simulation of an unbound wall in the x direction would be negligible. Our results showed that, at small Reynolds numbers, the simulations were insensitive to these lengths, so we typically chose 20D as the periodic repeat length. A particle was situated in the cell, as shown in Fig. 1, and a velocity distribution was assigned to the fluid boundary nodes representing a specific potential doublet in the particle.

Fig. 1. Schematic illustration of the flow field and the coordinate system.

Figure 2 shows the translation and rotation of a self-propelled particle in the central part of a fluid field. The streamlines depict the two typical swimming modes that occur under low-Reynolds-number flow.

Fig. 2. Illustration of streamlines for a single self-propulsive particle swimming in the channel: (a) translation, (b) rotation.
3.2. Translation of single particle in the channel

Compared with the boundary element method that Ishikawa et al. [19–22] proposed for use in simulating active suspensions, the method we propose here can reduce the huge amount of computational effort required as a result of the large number of singular numerical integrations involved. The advantage of LBM is that it avoids the need to solve Poisson’s equation and is easy to employ using parallel computation, which greatly reduces the computational resources required. On the other hand, unlike previously proposed computational schemes that lacked the ability to address complicated boundaries, the proposed method can construct the exact velocity distributions of particles with complex geometric shapes by arranging the basic singularities of different types of Stokes flow within a certain law. That is to say, the method we propose here possesses the advantages of both the LBM and the method of singularity distribution.

A self-propelled particle can swim continuously or at regular intervals. We call the interval between two adjacent times the switching time. We define the Reynolds number as a dimensionless constant associated with translational motion. Figure 3 shows the velocity distributions at different switching times. We can see that the distributions illustrate a dependence on the switching time. The different velocity distributions in the fluid depart from a continuous-swimming Maxwellian distribution in the areas of lower and higher velocities. A self-propelled particle results in a system that is out of equilibrium where energy is pumped into the fluid field, creating thicker high-velocity distribution tails. In addition, although the size of the particle is much smaller than the fluid domain, it can create thicker low-velocity distribution tails. The nearly identical low-velocity distributions with decreasing velocity are consistent with one of the Stokes flow characteristics, i.e., that the flow disturbance attenuation is slow and broadly distributed.

Fig. 3. (color online) Velocity distributions at different switching time. The different numbers stand for different interval time steps, for example, the number 100 means 100 (100×1 s). is the given translational velocity of the fluid boundary node where the outer normal direction of the particle is consistent with the particle movement direction, D is the particle length scale and υ is the kinematic viscosity. For Maxwellian, . v 0 is the fluid node velocities scaled by the average value of the fluid node velocities. is the normalised fluid nodes velocity distribution.

The velocity distributions for the translation of a single particle at different values are shown in Fig. 4, which is analogous to Fig. 3.

Fig. 4. (color online) Velocity distributions for different . For Maxwellian, .

The potential doublet at different positions along the orientation vector and the velocity distribution of the fluid nodes at different eccentric distances are shown in Fig. 5. We can see that changes in the eccentricity have no significant effect on the velocity distribution.

Fig. 5. (color online) Schematic illustration of a particle with eccentric potential doublet position and velocity distributions at different potential doublet eccentric distances. . For Maxwellian, .

The decrease in the velocity of a self-propelled particle when its activity is switched off was also studied at three Reynolds numbers. The results show that the shapes of the decay curves are consistent with a double exponential function. Figure 6 shows a fitted curve for the case of . This curve also coincides with all points at and . The regression line is of the following form:

(17)
with , , , , and .

Fig. 6. (color online) Particle translational velocity decay under different with time t after the driving is switched off. is the particle translational velocity, and is the particle translational velocity before the activity was switched off.

The timescales t 1 and t 2 stem from different dissipative mechanisms. The short timescale t 2 stems from the friction generated by relative motion between the particle boundary and the fluid, and the long timescale t 1 stems from the fluid drag force that exists before the particle becomes stationary.

Different eccentric distances can also lead to a various self-propulsion velocities, even though we assigned an equal-intensity potential doublet to the particle. A fitted parabolic curve of the following form is shown in Fig. 7 for different eccentric distances:

(18)
with , , , , , and .

Fig. 7. (color online) The change of velocity for the single particle with different , and is the average velocity of the self-propulsion particle.

As the figure shows, the eccentric distance can determine the self-propulsion efficiency. The farther the eccentric point is from the centre, the higher the particle velocity is.

3.3. Rotation of single particle in the channel

To explore the characteristics of the rotation of a single self-propelled particle in a channel, we defined a rotation Reynolds number and simulated rotational motion in the flow field by assigning a rotation speed to the fluid boundary of the particle. Three nearly overlapped distribution curves for three different values are shown in Fig. 8. For , there is a peak in the low-velocity region, and the curve is smooth in other regions. At very high velocities, the distribution decays so slowly as to be almost uniform when the self-propelled particle is rotating. The peak in the low-velocity region may be attributed to the effect of the boundary limitation which prevents the propagation and attenuation of fluid velocities toward further distances.

Fig. 8. (color online) Velocity distributions for different . is the given rotation velocity of the self-propulsion particle at fluid boundary. is the angular velocity of the particle, D is the particle length scale and υ is the kinematic viscosity. v 0 is the fluid node velocities scaled by the average value of the fluid node velocities.

The decay process of the rotation of a single particle is different from that of its translation, in that rotation requires less time than translation. Figure 9 shows a fitted curve for the case of . This curve also coincides with all points at and . A double exponential curve of the following form was used to describe the decay process and is shown in Fig. 9:

(19)
with , , , and .

Fig. 9. (color online) Particle rotation velocity decay under different with time t after the driving is switched off. is the particle rotational velocity, is the angular velocity of the particle, D is the particle length scale. is the particle rotational velocity before the activity was switched off.

The short timescale t 1 and the long timescale t 2 are determined by the friction and drag force, respectively.

3.4. Translation and rotation of single particle in the channel

A self-propelled micro-organism often switches its mode of motion between translation and rotation. For example, a bacterium cell usually exhibits a run-and-tumble motion in which it changes the direction of swimming by the rotation of any motor, i.e., unravelling of flagellar bundles. [15] Thus, it makes sense to study fluid velocity distributions after the sudden movement pattern changes of a single particle. Figure 10 shows that at , a particle abruptly changes its swimming direction by rotating anticlockwise to different degrees. The curves in the figure exhibit different peak velocities with increasing rotation angle. The greater the angle is, the smaller the peak velocity is, and the higher the distribution peak value is. The result is gradually upturned tails. At the beginning of an angle of , the effect of the translation velocity becomes insignificant because of the sufficiently long decay time. Rotational motion thus plays a major role in the distribution of velocity.

Fig. 10. (color online) Fluid node’s velocity distributions for different rotation angles at and .

A single particle rotating through a specific angle for various values was studied. The results are shown in Fig. 11. We can see that the velocity distributions in the high-velocity regions are coincident, but in the low-velocity region, the distributions are different. The curves shown in Figs. 10 and 11 resemble long-tailed non-Maxwellian distributions, which can be considered to be the result of the kinetic energy imparted to the fluid by the particle’s activity.

Fig. 11. (color online) Fluid node’s velocity distributions for different at .
4. Conclusion

A model that combines the lattice Boltzmann method with the singularity distribution method is proposed to simulate a self-propelled particle swimming (exhibiting translation and rotation) in a channel. Periodic boundary conditions were used in the simulations conducted, and the interaction between a single particle and the channel was not taken into account. The results show that the velocity distributions for translation and rotation of a single self-propelled particle deviate from a Maxwellian distribution, with high-velocity tails being observed because of the pumping of energy into the system. The effects of different potential doublet eccentricities on the translation velocity of a single self-propelled particle are significant. The double-exponential distributions decay when the swimming particle is switched off, which confirms that friction and drag force are the two major sources of resistance. The velocity distribution of a self-propelled particle is different from the Maxwellian distribution of a passive particle. No large differences in the velocity distribution were observed for different translation Reynolds numbers, rotation Reynolds numbers, or regular intervals.

Reference
1 Dreyfus R Baudry J Roper M L Stone H A Fermigier M Bibette J 2005 Nature 437 862
2 Yu T S Lauga E Hosoi A E 2006 Phys. Fluids 18 091701
3 Brennen C Winet H 1977 Annu. Rev. Fluid Mech. 9 339
4 Ajdari A Stone H A 1999 Phys. Fluids 11 1275
5 Paxton W F Kistler K C Olmeda C C Sen A Angelo S K S Cao Y Mallouk T E Lammert P E Crespi V H 2004 J. Am. Chem. Soc. 126 13424
6 Golestanian R Liverpool T B Ajdari A 2005 Phys. Rev. Lett. 94 220801
7 Gouin E Welch M D Cossart P 2005 Curr. Opin. Microbiol. 8 35
8 Leshansky A M 2006 Phys. Rev. E 74 012901
9 Schmittmann B Zia R K P 1995 Statistical Mechanics of Driven Diffusive Systems, in Phase Transitions and Critical Phenomena Vol. 17 New York Academic 1 p 1
10 Marchetti M C Joanny J F Ramaswamy S Liverpool T B Prost J Rao M Simha R S 2013 Rev. Mod. Phys. 85 1143
11 Wu X L Libchaber A 2000 Phys. Rev. Lett. 84 3017
12 Kim M J Breuer K S 2004 Phys. Fluids 16 L78
13 Hernández-Ortiz J P Stoltz C G Graham M D 2005 Phys. Rev. Lett. 95 204501
14 Underhill P T Hernández-Ortiz J P Graham M D 2008 Phys. Rev. Lett. 100 248101
15 Donald L K Ganesh S 2011 Annu. Rev. Fluid Mech. 43 637
16 Lighthill M J 1952 Comm. Pure Appl. Math. 5 1091
17 Blake J R 1971 J. Fluid Mech. 46 199
18 Magar V Goto T Pedley T J 2003 J. Mech. Appl. Maths. 56 65
19 Ishikawa T Simmonds M P Pedley T J 2006 J. Fluid Mech. 568 119
20 Ishikawa T Pedley T J 2007 J. Fluid Mech. 588 399
21 Ishikawa T Pedley T J 2007 J. Fluid Mech. 588 437
22 Evans A A Ishikawa T Yamaguchi T Lauga E 2011 Phys. Fluids 23 111702
23 Underhill P T Hernández-Ortiz J P Graham M D 2008 Phys. Rev. Lett. 100 248101
24 Saintillan D Shelley M J 2007 Phys. Rev. Lett. 99 058102
25 Zhang T Shi B C Chai Z H 2015 Acta Phys. Sin. 64 154701 (in Chinese)
26 Nie D M Lin J Z 2010 Chin. Phys. Lett. 27 104701
27 Li K Zhong C W 2015 Chin. Phys. B 24 050501
28 Song B W Ren F Hu H B Huang Q G 2015 Chin. Phys. B 24 014703
29 Zhang Y Pan G Huang Q G 2015 Acta Phys. Sin. 64 184702 (in Chinese)
30 Ladd A J C Colvin M E Frenkei D 1988 Phys. Rev. Lett. 60 975
31 Ladd A J C 1994 J. Fluid Mech. 271 285
32 Ladd A J C 1994 J. Fluid Mech. 271 311
33 Aidun C K Lu Y 1995 J. Stat. Phys. 81 49
34 Aidun C K Lu Y Ding E 1998 J. Fluid Mech. 373 287
35 Ding E Aidun C K 2003 J. Stat. Phys. 112 685
36 Ramachandran S Sunil Kumar P B Pagonabarraga I 2006 Eur. Phys. J. E 20 151
37 Alarcón F Pagonabarraga I 2013 J. Mol. Liq. 185 56
38 Chwang A T Wu T Y T 1974 J. Fluid Mech. 63 607
39 Chwang A T Wu T Y T 1975 J. Fluid Mech. 67 787
40 Chwang A T 1975 J. Fluid Mech. 72 17
41 Ku X K Lin J Z 2009 Phys. Scr. 80 025801
42 Bhatnagar P L Gross E P Krook M 1954 Phys. Rev. 94 511
43 Chen H Chen S Matthaeus W H 1992 Phys. Rev. A 45 R5339
44 Chen S Doolen G D 1998 Annu. Rev. Fluid Mech. 30 329
45 Happel J Brenner H 1965 Low Reynolds Number Hydrodynamics 1st edn Hague Martinus Nijhoff Publishers
46 Keller S R Wu T Y 1977 J. Fluid Mech. 80 259